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ABSTRACT 

We show that dark matter halos, in n-body simulations, have a boundary layer (BL), which neatly 
separates dynamically bound mass from unbound materials. We define T(r) and W(r) as the differential 
kinetic and potential energy of halos and evaluate them in spherical shells. We notice that, in simulated 
halos, such differential quantities fulfill the following properties: (i) The differential virial ratio 1Z = 
—2T/W has at least one persistent (resolution independent) minimum f, such that, close to f, (ii) the 
function w = —d\ogW/d\ogr has a maximum, while (iii) the relation 71(f) ~ w(f) holds. BL's are set 
where these three properties are fulfilled, in halos found in simulations of TCDM and ACDM models, 
run ad-hoc, using the ART and GADGET codes; their presence is confirmed in larger simulations of the 
same models, with a lower level of resolution. Here we find that ~ 97 % of the ~ 300 clusters (per model) 
we have with M > 4.2 • 10 14 fo — 1 Mq, owns a BL. Those clusters which appear not to have a BL are seen 
to be undergoing major merging processes and to grossly violate spherical symmetry. The radius f = r c 
has significant properties. First of all, the mass M c it encloses almost coincides with the mass M c i yn , 
evaluated from the velocities of all particles within r c , according to the virial theorem. Also, materials 
at r > r c are shown not be in virial equilibrium. Using r c we can then determine an individual density 
contrast A c for each virialized halo, that we compare with the "virial" density contrast A„ ~ 178 fi^' 45 
(where Q m is the matter density parameter) obtained assuming a spherically symmetric and unperturbed 
fluctuation growth. As expected, for each mass scale, A„ is within the range of A c 's. However the spread 
in A c is wide, while the average A c is ~ 25% smaller than the corresponding A„. We argue that the 
matching of properties derived under the assumption of spherical symmetry must be a consequence of 
an approximate sphericity, after violent relaxation destroyed features related to ellipsoidal non-linear 
growth. On the contrary, the spread of final A c is an imprint of the different initial 3-D geometries of 
fluctuations and of the variable environment during their collapse, as suggested by a comparison of our 
results with Sheth & Tormen analysis. 

Subject headings: galaxies: clusters: general — cosmology: theory — dark matter — large-scale 
structure of the Universe — methods: N-body simulations 



1. INTRODUCTION 

The inner region of dark matter halos in n-body sim- 
ulations is certainly virialized. A precise detection of the 
boundaries of such a region, however, is a long-lasting 
and not yet completely solved problem. If we try just to 
identify a "collapsed" region, i.e. a region detached from 
the cosmological expansion, the task is easier (see, e.g., 
Monaco & Murante 1999 and references therein). How- 
ever, in general, a collapsed region is not yet in virial equi- 
librium. The point is that, in a cosmological context, each 
structure is embedded in a more or less homogeneous back- 
ground, and we pass gradually from the inner halo to ex- 
ternal materials. The aim of this work is to show a precise 
criterion able to find virialized regions in dark matter ha- 
los. This criterion will be tested on simulations and found 
to be more efficient than previously adopted rules. 

As a matter of fact, several global properties of galaxy 
clusters, such as mass (M), radius (r) and, also, temper- 
ature (Tx) or luminosity (Lx) are strictly linked to the 
achievement of equilibrium. Analytical predictions on such 
global quantities were performed assuming that structures 
evolve from uniform, spherically symmetric fluctuations, 



preserve such symmetries during their growth and feel no 
influence from their environment (Gott & Rees 1975, Pee- 
bles 1980). However, even under such simplifying condi- 
tions, a cosmological constant is a significant complication. 
Lahav et al (1991) studied the growth of uniform spher- 
ical fluctuations within such context. They also outlined 
that the most direct way through which observations can 
distinguish among flat models with different Q m (matter 
density parameter) is comparing cluster abundances at low 
and high redshift. The dynamics of the spherically sym- 
metric growth of a uniform density fluctuation, in a model 
with Q m +Q\ = 1 (where is the vacuum density param- 
eter), was further explored by Eke et al (1996). Brian and 
Norman (1998) gave a polynomial expression for the ex- 
pected density contrast of virialized clusters, which reads 
A„ = A+Bx+Cx 2 , withx = fl m -l, and A = I8ir 2 ~ 178, 
B = 82, C = —39, and holds at any redshift z; this relation 
is well approximated by the expression: 

A„ = Afi°- 45 - (1-1) 

Under the same assumptions, using a Press & Schechter 
(1974) approximation, one can predict the expected mass 
function of clusters. 



department of Physics G. Occhialini, Universita di Milano-Bicocca, Piazza della Scienza 3, I 20126 Milano (Italy) 

2 I.N.F.N., Via Celoria 16, 120133 Milano (Italy) 

3 I.N.A.F., Osservatorio Astronomico di Torino - Torino (Italy) 



2 



Mass of clusters in simulations 



The impact of this approximation on the estimate of 
cluster global quantities could be appreciated only through 
numerical simulations. Much work was therefore devoted 
to this aim (see, e.g., Lacey & Cole, 1993,1994; Cole & 
Lacey 1996; Carlberg et al., 1996; Eke et al. 1996; Eke 
et al. 1998; Brian & Norman, 1998; Gardini et al. 1999). 
However, in the very analysis of the numerical work, the 
above approximations were often implicit. An example is 
the identification of dark halos based on spherical overden- 
sity (SO) algorithms, with reference to the density contrast 
A„ given in eq. (1.1). In this work we shall also use an SO 
algorithm (described below) as a selection criterion. How- 
ever, using eq. (1.1), to work out values for the cluster 
radius (r v ) or mass (M„), implies a bias. The values (1.1) 
are obtained under precise restrictions and numerical work 
should also verify their impact. In a sense, the problem is 
even more severe when the values (1.1) are used to define 
virializcd halos without requiring explicitly that they are 
spherical. 

In this work we plan to find a rule, suitable for defining 
virialized systems, taking into account that a variety of 
growth histories led to different final features of virializcd 
structures. A number of recent papers (Sheth & Tormcn 
(1999, 2002), Sheth, Mo & Tormen, 2001) had a similar 
motivation. They tried to study the growth of primeval 
fluctuations taking into account the effect of an ellipsoidal 
collapse. Then, in order to avoid having the collapse on 
some axis going to zero, they set up a recipe to freeze it 
out once it has shrunk by some critical factor. This freeze- 
out radius is explicitly chosen so that the density contrast 
at virialization for the whole halo is again A„ as given in 
eq. (1.1), i.e. in spherical growth models. In contrast, in 
this work the use of A„ to define virialized systems is com- 
pletely overcome. The point is that the detailed growth 
history of individual halos, in which different ellipticities 
and velocity anisotropies played an important role, is not 
expected to leave major traces on final halo shapes af- 
ter violent relaxation has occurred during the virialization 
process (Navarro, Frenk & White 1997). On the other 
hand, it may result in a significant spread of the equilib- 
rium density contrasts A c of individual halos. Therefore, 
when the virial radius of relaxed objects is sought, devi- 
ations from spherical symmetry and velocity anisotropies 
may not be as important as deviations of A c from A v . 

We expect that the value A„, given by eq. (1.1), is within 
the range of the actual A c values, but we think that it is 
important to first test how extended this range is, and 
how A„ is set with respect to actual values. Furthermore, 
a different A c implies a different radius and mass (r c and 
M c ) for each object. Accordingly, we expect that some 
individual object masses increase and others decrease. We 
will explore the extent of these changes on galaxy clusters 
in numerical simulations of critical CDM models with and 
without a cosmological constant. Such changes may have 
an impact on the mass function (see section 5 and Fig. 13 
and 14), while some scaling relations might be significantly 
altered, for instance the relation between M and Tx- 

In order to implement our program, we need a simple 
and effective rule to identify virialized material in clus- 
ters. In trying to find such rule we discovered a precise 
regularity in the transition region between material be- 
longing to clusters and the surrounding medium. In fact, 
the volume occupied by virialized materials is confined by 



a boundary layer (BL), whose depth Ar ranges around 50- 
100/i _1 kpc, and whose properties fulfill precise analytical 
prescriptions. Such a BL has been found in ~ 97 % of the 
clusters in our simulated models and is one of the findings 
of this work. 

The BL is not a physical confinement barrier, and par- 
ticles can travel outwards and inwards through it. On the 
contrary, it confines a volume, where material fulfills pre- 
cise equilibrium conditions, which cease to hold outside it. 
Inside the BL itself, the formal condition for virial equi- 
librium of an isolated system holds, as pressure forces at 
its upper and lower borders cancel each other; on the BL, 
kinetic energy has a minimum with respect to potential en- 
ergy, and matter density is rapidly decreasing. These con- 
ditions follow from a precise analytical requirement (Rw 
requirement) derived under the assumption of spherical 
symmetry. This Rw requirement will be only approxi- 
mately satisfied by halos in simulations, where the halo ge- 
ometries apparently violate spherical symmetry. We show, 
however, that normal deviations from a spherical shape in 
DM halos do not prevent them from having a BL. 

We explicitly note that we do not force the BL to exist, 
but we apply the Rw requirement and find it. For this 
reason, finding the BL is different from finding a given 
spherical overdensity; the latter always exists when aver- 
age densities of halos decrease with increasing radii, and 
this does not imply even an approximately spherical sym- 
metry. On the contrary, we could have bound and relaxed 
regions which are not confined by a BL: our analysis of nu- 
merical simulations shows that this is not the case. There 
is also a major difference between the Rw requirement we 
are defining here, and other prescriptions, such as SKID, 
which defines a halo as the set of all particles gravitation- 
ally bound to the halo itself. In fact, bound particles are 
not guaranteed to belong to virialized zones and may be- 
long to regions which are bound but still unrelaxed. 

In this paper we shall first discuss when the BL should 
exist and introduce the Rw requirement used to detect 
it. To have significant statistics based on this rule and 
to confirm its validity requires simulations performed with 
a large dynamic range, as can only be done by parallel 
computing. 

Finding and situating BL's allows a better evaluation 
of individual cluster masses, radii and density contrasts, 
as well as an analysis of their possible dependence on var- 
ious physical parameters. The study of the structure of 
dark matter halos in numerical simulations can also be 
put on a sounder footing. In previous analyses (see, e.g., 
Cole & Lacey 1996) the r dependence of the integral virial 
ratio — 2T(< r)/W(< r) was discussed and found not to 
be "useful for defining the boundaries of the virialized re- 
gion" . One of the findings of our analysis is precisely the 
possibility of detecting such a boundary, using the virial 
ratio for shells instead of spheres. Fig. 1 compares integral 
and differential virial ratios, showing how much more infor- 
mation appears to be contained in the latter (see Sec. 3 for 
details on how lZ(r) points are calculated) . Previous anal- 
yses have also stressed halo properties apparently depend- 
ing "on how groups are identified" (Sheth et al., 2001). 
The existence of a border, with precise physical proper- 
ties, secluding cluster material lets one overcome any such 
ambiguity. 

It may be also important to remember that according 
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to the definitions given above, in this work A„, r v , M v 
are quantities worked out starting from the fixed density 
contrast given by (1.1), while A c , r c , M c are the same 
quantities when worked out from the setting of the BL 
around each cluster. Let us also point out that CHb the 
procedure discussed in this work is not the basis for an 
CHe alternative cluster-finding algorithm, such as, e.g., 
Friends-of- Friends, SKID or SO. On the contrary, our anal- 
ysis assumes that cluster locations are given. Here we give 
a different prescription for defining the physical dimension 
of the cluster, not its position. 

The plan of the paper is as follows. In section 2 we give 
suitable information on the simulations used, which were 
partially run ad-hoc, for this work. In section 3 and 4 
we define the boundary layer (BL) and discuss how it can 
be found in simulated clusters. In section 5 we show the 
results of our work, and in section 6 we discuss them and 
some future perspectives. 

2. MODELS AND N-BODY SIMULATIONS 

We use two main sets of simulations, performed with dif- 
ferent codes: the parallel AP3M N-body code described 
in Gardini et al. (1999), the parallel PM ART code de- 
veloped by Kratshov & Klypin (1997), and the parallel 
tree-code GADGET (Springel et al. 2001). 

The code of Gardini et al. (1999) was developed from 
the serial public AP3M code of Couchman (1991), extend- 
ing it to different cosmological models and allowing suit- 
able flexibility in particle masses. Using this code, two sim- 
ulations were performed. The first of them, dealing with 
a "tilted" Einstein-de Sitter model (hereafter TCDM), is 
the same already considered by Gardini et al. (1999); the 
normalization of the run was, however, rescaled to yield 
£78 = 0.55 at the final step, instead of erg = 0.61, so that 
the final abundance of rich clusters is the same for both 
models at the final epoch (as ususal, erg is the m.s. density 
fluctuation on the scale of 8 /i _1 Mpc, h being the Hubble 
parameter in units of 100 km/s/Mpc). The latter simu- 
lation dealt with a ACDM model, i.e. a flat CDM uni- 
verse with non-zero cosmological constant with o~g = 1.08. 
Owing to recent observational results (see, e.g., Schuecker 
et al, 2002), these normalizations are both rather high, 
as a reasonable value for erg in ACDM models is <~ 0.75; 
however, this has no impact when studying inner cluster 
dynamical properties. 

The above simulations (denoted A and B, respec- 
tively) are performed in cubes with sides of 360 h~ 1 Mpc. 
CDM+baryons are represented by 180 3 particles, whose 
individual mass is 2.22 • 10 12 /i _1 Af© for TCDM and 0.777- 
1O 12 /i _1 M for ACDM. We use a 256 3 grid to compute 
the FFT's needed to evaluate the long range contribution 
to the force (PM) and we allow for mesh refinement where 
the particle density attains or exceeds <~ 30 times the mean 
value. The starting redshifts are = 10 for TCDM and 
Zi n — 20 for ACDM. The particle sampling of the density 
field is obtained by applying the Zel'dovich approximation 
(Zel'dovich 1970; Doroshkevich et al. 1980) starting from 
a regular grid. We adopt the same random phases in both 
A and B. 

The number of steps were 1000 equal p-time steps (the 
time parameter is p oc a 2 / 3 , where a is the expansion fac- 
tor). The comoving force resolution, given by the softening 



length, is r\ ~ 112/i kpc, yielding a Plummer equivalent 
softening e p i ~ 40.6 /i kpc. These simulations were de- 
scribed in more detail in Gardini et al. (1999), and were 
also used in Maccio et al. (2002). The parameters of the 
models are reported in detail in Table 1. 

In order to see how resolution affects the results of 
the Rw requirement, we performed further simulations 
with the ART code (courtesy of A. Klypin). This al- 
lows us to select regions inhabited by clusters and to re- 
run them with increasing particle mass resolution: from 
7.7 x lO 11 fi m /i- 1 M to 1.2 x lO^On^MQ. 

The ART code (Adaptive Refinement Tree: Kravtsov et 
al. 1997, Knebe et al. 2000) starts with a uniform grid, 
which covers the whole computational box. This grid de- 
fines the lowest (zeroth) level of resolution of the simu- 
lation. The standard Particles-Mesh algorithms are used 
to compute the density and gravitational potential on the 
zeroth-level mesh. The code then reaches high force res- 
olution by refining all high density regions using an auto- 
mated refinement algorithm. The refinements are recur- 
sive: the refined regions can also be refined, each subse- 
quent refinement having half of the previous level's cell 
size. This creates a hierarchy of refinement meshes of dif- 
ferent resolution, size, and geometry covering regions of 
interest. Because each individual cubic cell can be refined, 
the shape of the refinement mesh can be arbitrary and 
effectively match the geometry of the region of interest. 
This algorithm is well suited for simulations of a selected 
region within a large computational box. 

The criterion for refinement is the local density of parti- 
cles: if the number of particles in a mesh cell (as estimated 
by the Cloud-In-Cell method) exceeds the level n t hresh, the 
cell is split ("refined") into 8 cells of the next refinement 
level. The refinement threshold may depend on the refine- 
ment level. The code uses the expansion parameter a as 
the time variable. Besides spatial refinements, during the 
integration time refinements are also performed. 

The ART code can handle particles of different masses. 
This lets us increase the mass (and correspondingly the 
force) resolution of regions centered around the highest 
mass clusters. The multiple mass resolution is imple- 
mented in the following way. We first set up a realiza- 
tion of the initial spectrum of perturbations, so that ini- 
tial conditions for our largest number (512 3 ) of particles 
can be generated in the simulation box. Coordinates and 
velocities of all the particles are then calculated using all 
waves ranging from the fundamental mode k = 2tt/L to 
the Nyquist frequency k = 2n/L x N 1 ^ 3 /2, where L is 
the box size and N is the number of particles in the sim- 
ulation. We obtained lower resolution regions by merg- 
ing high resolution particles into particles of larger mass 
where needed. This process can be repeated to get still 
lower resolution. The larger mass (merged) particles have 
velocities and displacements equal to the average of the 
velocities and displacements of the smaller-mass particles. 
Using this technique we firstly generated initial conditions 
for low-resolution simulations. These simulations were run 
from their initial redshift to z = and they were used to 
locate the Lagrangian regions forming the most massive 
clusters. Then, we were able to obtain initial conditions 
at high resolution (both in force and mass) in the selected 
Lagrangian zones, with the same random Fourier phases. 
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Such zones were surrounded by lower and lower resolu- 
tion shells taking into account the tidal field which acts 
on them. Finally, these ICs were used to run the high- 
resolution cluster simulations. 

The simulations employed here were performed us- 
ing 128 3 zcroth-level grid in a computational box of 
180/i _1 Mpc. The threshold for cell refinement (sec above) 
was low on the zeroth level: nthrosh(O) = 2. Thus, ev- 
ery zeroth-level cell containing two or more particles was 
refined. The threshold was higher on deeper levels of re- 
finement: n t hrosh = 3 and n t hrosh = 4 for the first level and 
higher levels, respectively. 

For the low resolution runs the step in the expansion 
parameter was chosen to be Aao = 2 x 1CP 3 on the zeroth 
level of resolution. This gives about 500 steps for particles 
located in the zeroth level for an entire run to z = and 
128.000 for particles at the highest level of resolution. 

Using the ART code we performed two simulations: a 
ACDM model and a TCDM model, that we will indicate as 
C and D, respectively. These models are similar, although 
not identical, to those of Gardini et al. 1999; their parame- 
ters are listed in Table 2. We identified the clusters in these 
zero-level simulations and selected the 6 largest mass clus- 
ters for each cosmological model. They were re-run with 
increased mass (and force) resolution, using the method 
described above. For every selected cluster we have three 
runs with resolutions 4 of 128 3 (particles mass M p = 
7.7 x 10 n f7 m h- l M Q ), 256 3 (M p = 9.6 x 10 10 O m /i" 1 M ) 
and 512 3 (M p = 1.2 x 10 10 f7 m /i" 1 M Q ) with more than 
350.000 particles within a radius of 2.0 Mpc h . 

We then run 3 out of the 6 high resolution ACDM clus- 
ters with the public parallel tree-code GADGET (Springcl 
et al. 2001), to verify that our results do not depend upon 
the peculiarities of one N-body code. One of the clus- 
ters was magnified with both ART and GADGET. The 
initial conditions were set as described above. We used 
a Plummer-equivalent softening length e = 10.98; this is 
also the linear dimension or our highest-resolution ART 
cell. We chose the time-step criterion 3, based on the local 
dynamical time, with a tolerance of the integration error 
Eint = 0.2, which gave us more than 100000 time steps 
from Zi n to z = 0. The criterion for opening the tree-cells 
is based on the absolute truncation error in the multipole 
expansion, with a tolerance on the error on the force of 
Ep = 0.02 (see Springel et al. 2001 and the code user 
manual for more details). 

In all of our numerical simulations, clusters were identi- 
fied using a SO algorithm. As a first step, candidate clus- 
ters are located by a FoF algorithm, with linking length 
A = (f>x d (here d is the average particle separation), keep- 
ing groups with more than Nf particles. We then perform 
two operations: (i) we find the point, Cw, where the grav- 
itational potential, due to the group of particle, is mini- 
mum; (ii) we determine the radius, f, of a sphere centered 
in Cw where the density contrast is A„. Using all parti- 
cles in this sphere we perform again the operations (i) and 
(ii). The procedure is iterated until we converge onto a 
stable particle set. The set is discarded if, at some stage, 
we have less than Nf particles. If a particle is a potential 
member of two groups it is assigned to the more massive 
one. (Gardini et al. 1999 also describe this SO algorithm 



in more detail and compare it with other group identifica- 
tion algorithms, e.g. Governato et al. 1999). In this work 
we set 4> = 0.2 and take Nf so to have a mass threshold 
3.0 x 10 13 /i _1 M Q . Above this mass threshold there are 
~ 10000 halos in the simulations A,B and ~ 1250 halos in 
the simulations C,D. 

Tests performed by Gardini et al. (1999) show that no 
cluster is missed above 4.2 x 10 14 ft, _1 M Q . We have 303 
such clusters in A, 316 in B, 40 in C, 36 in D; 6 clusters of 
C and D where then blown up to various resolution levels. 



3. 



THE VIRIAL THEOREM APPLIED TO SINGLE 
CLUSTERS: THEORY 



Finding the volume where cluster materials are in virial 
equilibrium may seem a tough and somewhat ambiguous 
task. If kinetic and potential energies, within a radius r, 
are defined according to: 



2T(<r)=£ 



(r«<r) 



W(< r) 



Cm 2 



;<r) 



the virial ratio 



TZ(< r) = 



2T(< r) 
~W{< r) 



(3.1) 



(3.2) 



(3-3) 



should be unity, in a virializcd sphere of radius r, pro- 
vided that it is fully isolated. On the contrary, in the real 
world, as well as in simulations, the virialized volume is 
bordered by infalling and outgoing material, possibly trav- 
eling through a partially depleted area. For this reason, 
the convergence of 1Z(< r) to unity, at a precise radius r, 
should not even be expected. 

In principle, one can take into account such perturba- 
tions by fixing a radius r well inside the cluster volume 
and adding a pressure term to take into account external 
material, as has been done in observational work aiming 
to estimate the mass of optical clusters (see, e.g., Girardi 
et al., 1998). Accordingly, the virial equilibrium condition 
reads 

2T + W = 3Vp (3.4) 

(V is the volume occupied by virialized material and p 
is the pressure at its borders). Eq. (3.4) translates the 
problem into finding a reasonable estimate for p. 

The problem can also be considered from a slightly dif- 
ferent point of view: if we find, around the cluster, a bound- 
ary layer (BL), which is in virial equilibrium by itself, and 
above which virial equilibrium is lost, its very material and 
all material it borders, down to the cluster center, need to 
be in stationary equilibrium. Particles may pass through 
this BL in both directions, but the overall virial balance 
shall keep constant in time. 

We will now assume that this BL is spherical. The con- 
ditions found under this restriction are, however, well ap- 
proximated by simulated halos. This seems to indicate 
that, even after a strongly asymmetric growth, the final 

4 Since our averaging procedure corresponds to using different mesh sizes when generating ICs for a simulation, we will refer to the size of 
the higher resolution mesh used in each simulation IC when quoting resolutions. 
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virialization process tends to partially suppress elliptical 
features. 

The r dependence of the potential energy can be de- 
rived by performing the sums in cq. (3.2) starting with 
the outside: 



W(< = -E E 



Gm 2 



(3.5) 



1 (r*<r) 3 (r 3<ri ) 

= -E > 

* (n<r) 

where the last term is the very definition of i?. In the 
continuous limit (and still assuming spherical symmetry) 



Z(r) = Gm I" d 3 r' P ^ 



\r — r'\ 



(3.6) 



with a suitable density, p. Quite in general, a volume in- 
tegral of p(r) increases with r. Henceforth, 



Z(r) = C(r/f)- w , 



(3.7) 



where C is a normalization constant evaluated at an arbi- 
trary position f and w(< 1) will, however, depend on r. 
Accordingly, it must be that 

d? r 
r=—(r) = -\w + rw' In (-)]Z(r) . (3.8) 
dr r 

Let us now suppose that there exists an interval Ar = 
r + — r_ , which satisfies the following conditions: 

(i) it is in virial equilibrium, 

(ii) the r dependence of pressure is such that the r.h.s. 
of eq. (3.4), yielding a term of the form 3(p + V + — 
P-V-), can be neglected in the virial balance, 

(iii) inside Ar, w is constant. 

Further comments on the (ii) requirement can be found in 
Appendix A, where we show that halo profiles are expected 
to approach it. Together with (i), it implies that 



E mv 2 -r~(r t )=Q 

» (r.GAr) 



(3.9) 



in the layer. 

Let us then take into account the condition that w is 
constant and use eq. (3.8) with w' = and the virial ratio 



iz. 



(3.10) 



to obtain that 

1Z = w, (3.11) 

all along the interval Ar. This equation coincides with 
cq. (3.29) in Goldstein, Pole & Safko (2002) (in turn, 
the latter equation yields eq. (3.9) hereabove, in the limit 
w = const.). 

Requiring that w' = 0, then, also yields that 



dr 







(3.12) 



in Ar. Vice versa, if the eqs. (3.11)-(3.12) are both ful- 
filled in a layer of depth Ar, this layer is at rest and in 
virial equilibrium. We thus define boundary layer (BL) a 
region of depth Ar satisfying eqs. (3.11)-(3.12). 

Let us stress that the condition (iii) is essential to set 
the Rw requirement, defined by the eqs. (3.11)-(3.12). Un- 
less w is constant along a suitable interval, we can neither 
replace ridZ(ri)j dr by —wZ(ri) in eq. (3.9), nor factor- 
ize w out of the sum, so to obtain eq. (3.11). Moreover, 
should eq. (3.11) hold on a single r, we could have w' = 
there, with TZ' ^ 0. Finding a layer fulfilling the Rw re- 
quirement, goes therefore beyond finding a layer in virial 
equilibrium. As we shall see, however, once the Rw re- 
quirement is fulfilled, we have found a layer bounding a 
virialized region. 

Let us also stress that the Rw requirement does not 
imply that particles cannot pass through the BL; it only 
prescribes that such (possible) passages do not violate its 
stationary conditions. However, if we find about a cluster 
a boundary layer, it will act as a confinement barrier to 
inner kinetic and potential energies and mass. Notice that, 
in principle, inner materials might not be in virial equilib- 
rium themselves; even in this case, however, there can be 
no net exchange of mass between inside and outside and 
the possible energy exchange is limited to the work done 
by tidal torques due to anisotropies. 

Let us now show that no further material, outside the 
BL, can be in virial equilibrium. We will show that this 
is certainly forbidden if the values, taken by TZ and w at 
r + , are a minimum and a maximum, respectively, when 
compared with values at r > r + . This simple condition, 
however, is not strictly necessary; it is sufficent that the 
difference TZ — w increases, becoming positive, at r > r+. 
This will be shown here below. Let us however antici- 
pate that, while the Rw requirement can be fairly easily 
tested, it is not easy to deal with such further require- 
ment. The tests performed on a large number of clusters 
and described in the next sections, however, let us conjec- 
ture that, when the Rw requirement is fulfilled, this fur- 
ther requirement is somehow inescapable. We now show 
why, at r > r+, the equilibrium conditions are violated if 
a boundary layer is found. Notice that, at r > r + , we can 
still assume that W = W(r + )(r/r + )~ w , but w, starting 
from the value it had in Ar, will now vary with r and, 
according to eq. (3.8), 



TZ — w = w'rln(r/r + ) ; 



(3.13) 



this equation cannot be fulfilled, if the difference 1Z — w, 
vanishing up to r + , increases at r > r+, becoming posi- 
tive. In fact we required that w is maximum and therefore 
w' < 0. 

Let us finally notice that if a layer is characterized by 
higher w values, this means that it is emptier than inner 
or outer layers. A maximum of w, therefore, indicates a 
low density layer. In turn, a minimum of 1Z indicates a 
low kinetic energy, with respect to potential energy. BLs, 
therefore, are partially depleted r-intervals, where parti- 
cles, on average, are particularly slow. 
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4. THE VIRIAL THEOREM APPLIED TO SINGLE 
CLUSTERS: PRACTICAL USE 

The mathematical definition must now be applied to 
numerical simulations. A failure to find BLs could be at- 
tributed to insufficient numerical resolution or to bad vio- 
lations of spherical symmetry. However, in principle, even 
if the resolution is enough and no substantial departures 
from sphericity are visible, BLs could still be absent. 

Let us now describe our procedure. After finding clus- 
ters in simulations and their centers Cw, as described 
in sec. 2, we evaluated W(r) and TZ(r), for spherical r- 
intervals containing a fixed number 

A = M 4 /(2 4 m p ), (4.1) 

of particles (here M4 of the total mass within 4.5 /i _1 Mpc 
from Cw and m p is the particle mass). However, succes- 
sive points along r are obtained by shifting outwards by 
N/8 particles. In the relevant r range, the above criterion 
yields N ~ 10 2 for TCDM (simulation A) and N ~ 3 • 10 2 
for ACDM (simulation B), while successive points lay at a 
distance ~ 20 /i _1 kpc. This is below (but not much below) 
the force resolution (~ 100/i _1 kpc), as desired. 

We do not expect BL's to be much thicker than our force 
resolution; accordingly, we first seek the minima of TZ. Ow- 
ing to the procedure by which TZ points are worked out, 
a minimum is considered significant only if it is such with 
respect to 8 neighbors. Finding the maxima of w is harder, 
as w(r), obtained by differentiating W(r), should then be 
further differentiated to find its maxima. In practice, how- 
ever, this second order differentiation is not needed. Quite 
in general, there will be just a few points r m where TZ has 
significant minima. Let then lZ m be the virial ratio at such 
points and let us fit W(r) with the expression 

W = A(r/r m )- n ™ (4.2) 

in v = 8 points > r m , covering an interval Ar ~ 100- 
150/i _1 kpc. The fit is performed in two steps: we first 
minimize the function <p = W/W — 1 just using A as a 
parameter; then, in cq. (4.2), we allow lZ m to take val- 
ues different from the value of 1Z in r m , and perform a 
2-parameter fit. 

We considered the fit to be satisfactory, if two conditions 
hold: (i) When performing the first fit, the residual 

V 

X 2 = 5> 2 (r.) < Hr 3 ^ , (4.2) 

(i.e., a mean square discrepancy <~ 3%, between the nu- 
merical values of W(r) and its fitting expression, was our 
limit of tolerance). This detects r-intervals where eqs. 
(3.11)-(3.12) arc both fulfilled, (ii) When performing the 
second fit, the best-fit value of 7Z m is within 1—cr from the 
actual 7Z m value. Here we test that also the slope of W is 
not too far from what is required. Such constraints should 
therefore select points where TZ is minimum and (nearly) 
intersects w. 

As outlined above, one might reasonably expect that a 
large fraction of clusters, at their boundaries, is far enough 
from sphericity that the expected system of maxima and 
minima, as well as their coincidences, is disrupted. The nu- 
merical noise in calculating w is also quite large and can be 



expected to create severe problems. On the contrary, we 
found that all clusters with mass M >~ 4.2 • 10 14 /i _1 M Q 
have at least one value of r m fulfilling the above require- 
ments, while some clusters have 2 or 3 points where the 
fitting criteria are fulfilled. In the few cases, when the 
fitting criteria were satisfied at more than one r m , we dis- 
carded r m values yielding A c outside an interval 0.3 A„- 
3A„. This selection criterion left us with ~ 97% of the 
initial clusters. When more than one r m was still available, 
we selected the one corresponding to the smallest \ 2 ', this 
last criterion was needed for ~ 10% of the clusters only. 

Let us recall again that the requirement of intersection 
between a minimum of TZ and a maximum of w was found 
for spherical systems. When we deal with the halos of our 
simulations, sphericity is likely to be an approximation. 
When a single intersection exists, however, it is natural 
to guess that this is the point where the Rw requirement 
is best approached. However, both in this case and when 
several r m points are found, it is natural to perform a 
closer inspection to make sure that no unexpected feature 
may bias our understanding of the physical situation. 

This closer inspection requires an increased resolution. 
As described in sec. 2, a TCDM model and a ACDM model 
were then reconsidered within a smaller box with sides of 
180 h^Mpc, making use of the ART and GADGET codes. 
This gave us ~ 1/8 of the clusters provided by the AP3M 
simulations, but enabled us to focus on particular clusters, 
making use of the ART package facilities to increase the 
resolution there while still keeping identical boundary con- 
ditions. In order to test that the regularities found are not 
linked to peculiarities of the ART package, some clusters 
were magnified using GADGET, instead of ART. Alto- 
gether we have 12 clusters magnified, 6 for each model. 4 
of the latter were obtained using ART and 3 using GAD- 
GET; hence, one cluster was magnified with both codes 
showing that our results are independent of the n-body 
code used. In fact, the final positions of particles have 
just minor differences, while TZ and W coincide. A num- 
ber of subsequent tests also concerned the stability of the 
minima of TZ. In general, the number of local minima r m 
decreases when the resolution is increased, but the point 
where W fits the expression (4.1) with the worst resolution 
is a minimum also with greater resolutions. In just a few 
cases, new minima arise when the resolution is improved; 
but in no case is a minimum yielding the BL erased when 
changing resolution. The "worst" cases for each model are 
shown in Fig. 2. Here we have a TCDM cluster, taken 
from the simulation C, which passes from 3 minima (128 3 
particles), to 1 minimum (256 3 particles) and, again, to 3 
minima (512 3 particles). The BL is displaced outwards by 
~ 150 ft kpc, when we pass from 128 to 256; then, when 
passing from 256 to 512, no appreciable shift occurs. Out 
of the 12 clusters considered, only this one presents a shift 
slightly above the resolution of the initial simulation. The 
ACDM cluster of simulation D shown in Fig. 2, shows 4 
minima with 128 3 particles, 2 minima with 256 3 particles 
and, again, to 3 minima with 512 3 particles; but when 
changing resolution, the shifts are < 40/i~ 1 kpc. 

For these 12 clusters we follow in detail the joint be- 
havior of TZ and w. We show them in Fig. 3 for the best 
and worst cases we found. Of course, a precise coincidence 
between a minimum of TZ and a maximum of w would indi- 
cate that, within the resolution limits, the system is spher- 
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ical. Seeking this kind of coincidence, therefore is somehow 
fatuous. Furthermore, in order to work out w, we ought 
to differentiate W, which is obtained from a discrete point 
distribution. Hence, although the nice coincidence of the 
best case is quite appealing, this case is not safer than the 
worst case we show. By tolerating a 3 % disagreement in 
W, for r/r m <^ 1.1, we allow for a tolerance <~ 0.3 on 
the exponent w. In the worst case, we find a w-TZ discrep- 
ancy of this order, but in most cases, |w-7\L| turns out to 
be ~ 0.1. We tentatively conclude that this resolution is 
sufficient to show physical features and that the level of 
the w-TZ (dis)agreement is a measure of the a-sphericity 
of actual systems. The important issue is that we always 
find a maximum of w to match the minimum of TZ, quite 
close to the r m selected at the initial resolution level. 

This point can be further illustrated by Fig. 4, which 
shows a plot similar to Fig. 3, for a typical case at the 
initial resolution level. Both the TZ and the w curves are 
much noisier here. TZ, in particular, shows a wide set of lo- 
cal minima; the only persistent minima, however, are those 
marked by thick filled circles. As expected, the behavior 
of w is still noisier and most of the oscillations shown are 
clearly numerical. The role of w, however, just amounts to 
choosing among the points marked by a thick filled circle 
and selects the point marked by a vertical arrow. Even 
at this resolution level, the presence of a maximum of w, 
close to such point, can be suspected. When the resolution 
is increased, this maximum is more clearly visible, while 
no w maxima exist close to the other TZ minima. 

A further preliminary test of our technique, concerning 
galactic-sized and unvirialized DM halos, is given in Ap- 
pendix B. 

Before concluding this section, let us draw the reader's 
attention on the setting of the BL around particle agglom- 
eration, as shown in Fig. 5 for a ACDM cluster. This fig- 
ure shows that the BL is actually set outside of the main 
matter agglomeration; its spherical shape is somehow in 
contrast with the apparent a-sphericity of the cluster ma- 
terials and easily justifies minor discrepancies from results 
based on an assumption of full spherical symmetry. 

5. RESULTS 

Once the sphere confining cluster materials is set, we 
can evaluate the density contrast A c and the mass M c in- 
side it. In Figs. 6 and 7, points give A c and M c for all 
clusters in simulations A and B, respectively. They show 
a rather wide spread of A c values. By subdividing the 
M c abscissa in intervals of constant logarithmic width, we 
evaluate the average density contrast in each of them, to 
seek systematic trends with mass. Averages are still sub- 
ject to a significant uncertainty owing to the spread of A c . 
They are shown, at the 1-tr level, in the plots. At high 
masses, some logarithmic intervals of Fig. 7 are empty or 
contain a single object. In the latter cases no error bar is 
given. 

There seems to be no evidence of any peculiar trend of 
density contrasts with mass apart from, perhaps, a mod- 
est indication of an increasing density contrast at very 
high scales in ACDM. It is therefore possible to consider 
the overall average among A c 's. This average is indicated 
by the continuous horizontal line and compared with the 
"virial" density contrast A„, as given by eq. (1.1). In both 
cases, A„ (dotted line) is well inside the range of the den- 



sity contrasts we found; the average A c , however, in both 
cases, is smaller than A„ by ~ 25 %. 

It is, however, clear that gravitationally evolved, sta- 
tionary objects are not uniquely identified with a given, 
fixed density contrast. In further work, we plan to deal 
with the impact of variable A c on the physical character- 
istics of DM halos, such as concentration, average angular 
momentum, density and velocity profiles. It may also be 
significant to test our results against the effects of different 
particle selection criteria, e.g. using the SKID algorithm 
instead of the SO one to pre-select the DM halos. 

In Figs. 8 and 9 we plot r c against M c , for the simu- 
lations A and B, respectively. Here, error bars account 
for the variance of r c around its average, which is wider 
than the error of the average r c , shown in Fig. 6. Contin- 
uous and dotted curves show the expected behavior of r c 
against M c , when the density contrast is set either at the 
average A c value or at A„. 

It is also significant to compare the plot of r c against M c 
with the behaviour of the radii r c against the masses M v 
of the corresponding clusters. This comparison is shown 
in Fig. 10, for the simulation A only; in order to avoid 
confusion, we omit error bars which, however, are approx- 
imately of the same size in both cases. The highest mass 
points, for which error bars cannot be evaluated, are also 
omitted. Filled circles yield r c vs. M c , empty circles yield 
r c vs. M v . Filled circles nicely fit the line worked out 
from the average density contrast A c ; only at high masses, 
where the statistics is poorer, some of the circles are far 
from this line. The same line is also a reasonable fit for 
the open circles. These, however, are systematically far- 
ther at all mass scales. In the figure, however, we also 
plot an horizontal line showing the average of the r c val- 
ues of all clusters with mass > 4.2 • 10 14 /i _1 M Q . Even this 
average can be considered a reasonable fit for open circles 
(an unweighted % 2 evaluation gives similar probabilities to 
both fitting lines). This plot shows that there are actual 
dangers, if quantities worked out assuming a constant A„ 
are compared with real data. In this case, one might ten- 
tatively suggest that there is a typical, mass independent, 
cluster radius R ~ 1.56/i _1 Mpc, quite close to the Abcll 
value. Our procedure, instead, outlines that this arises be- 
cause some clusters were attributed a biased mass and the 
corresponding points were then set at a wrong abscissa. 
Although, altogether, this abscissa reshuffling is substan- 
tially casual, the danger that an characteristic scale erro- 
neously appears is far from absent. 

A critical result of our analysis is however shown in 
Figs. 11 and 12 (for simulations A and B, respectively). 
In simulations, the "numerical" cluster masses can be eas- 
ily evaluated by summing up particle masses. The actual 
approach, in the real world, ought to be quite different 
and is based, first of all, on an analysis of galaxy velocities 
or X-ray emitting gas temperature. These analysis make 
use of the virial hypothesis; therefore, measured masses 
of galaxy clusters are estimates of their dynamical masses. 
It is therefore important to test how "numerical" cluster 
masses in simulations agree with the masses 



evaluated averaging over the velocities of the particles 
within r„ ;C . This comparison has been often performed, 
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showing a reasonable coincidence between M v and Md vn . 
The results of such fit are shown in Figs. 11 and 12 by 
the dashed histograms, which confirm the slight excess of 
Mdyn vs. M v already noticed by previous authors. In the 
same plots we also show the results of a fit between Md yn 
and M c , obtained on the basis of the setting of the BL. 
There seems to be little doubt that the latter fit is better. 
The average values of Md yn /M c are both ~ 0.97 ± 0.03, 
quite close to unity and consistent with it at the \-a level. 

We have also compared the shapes of the cluster (inte- 
gral) mass functions n c (> M c ) and n c (> M v ). In Fig. 13 
and 14, upper and lower plots show their behaviors, respec- 
tively. In the same plots, we also show Press & Schechter 
and Shcth & Tormen theoretical curves. They can be ob- 
tained from two different expressions of 

M 

f{v)vd\ogv = —n c (M)M dXogM ; (5.2) 

Pm 

here p m is matter density, the bias factor v = 8 c /ctm is 
the ratio of the critical overdensity in the spherical growth 
model to the rms density fluctuation on the length scale 
corresponding to the mass M, and n c (M) is the differen- 
tial mass function, which shall then be integrated to obtain 
n c (> M). 

In the Press & Schechter formulation, assuming spheri- 
cal unperturbed fluctuation growth, 

f{v) v = v/27^ v exp(-i/ 2 /2) , (5.3) 

while Sheth & Tormen, taking into account the ellipsoidal 
nature of the collapse, argue that 

f(y) v = A(l + v' - 2q )^2/^ v exp(-i/ 2 /2) , (5.4) 

with v' = \fav\ this expression introduces the parameters 
A, q and a, however constrained by the requirement that 
all matter is bound in objects of some mass M, however 
small, and therefore the integral of f(y) must be unity. In 
our plots, A is normalized to yield the number of clusters 
found in simulations, and this constraint is automatically 
satisfied. Sheth & Tormen (1999) showed that, in the GIF 
simulations for an SCDM, OCDM and ACDM models (see 
Kauffmann ct al 1999), a good fit to the numerical mass 
function was obtained if q = 0.3, a = 0.707. 

In Fig. 13 and 14 we report the mass function obtained 
from eq. (5.4) for these values, for TCDM and ACDM, re- 
spectively. The fit with the n c (> M v ) is surely better than 
using eq. (5.3). We also tried to modify the parameter a, 
that, according to Sheth & Tormen arguments, is related 
to the density contrast of virialized structures. As a mat- 
ter of fact n c (> M c ) is derived with a variable density 
contrast, as explained in previous sections, and, however, 
with an average final density contrast which is <~ 25 % be- 
low A„. It is therefore reasonable to expect that a lower 
value of a is needed to compare theoretical predictions 
with the numerical behavior of n c (> M c ). In Fig. 13 and 
14, we show the expression (5.4), both for a — 0.707 and 
0.5. It is also evident that, while a = 0.707 provides a 
good fit for M„-mass functions, the fit obtained provided 
by a = 0.5 for the M c -mass function is perhaps better. 

The results of these comparisons, altogether, are that: 
(i) No major variations in the numerical mass function 



arises when M v is replaced by M c . (ii) Variations how- 
ever exist, which are of the same order of the differences 
between the standard theoretical expression (5.3; Press & 
Schechter) and the improved theoretical expression (5.4; 
Sheth & Tormen). (iii) When the latter expression is used, 
the fits to the numerical mass functions are best if differ- 
ent values of the parameter a are used for M v and M c 
(<~0.7 and 0.5, respectively). The differences between the 
mass functions n c (> M v ) and n c (> M c ) are more evident 
in the intermediate mass range (7 • 10 14 /i _1 M Q < M < 
10 15 /i _1 Mq). They arise from a widespread "re-shuffling" 
among the masses of DM halos, when the two definitions 
are used, causing thereby different numbers of halos in dif- 
ferent mass intervals. This can also be seen from Fig. 10, 
where the mean cluster radii (r c , r v ) of halos in this mass 
range show appreciable differences. 

6. CONCLUSIONS 

Observational work on galaxy clusters made a giant leap 
forward thanks to X-ray observations and, in turn, cluster 
data are providing more and more cosmological informa- 
tion. The ROSAT All Sky Survey (RASS) allowed the ex- 
traction of large samples of galaxy clusters, like REFLEX 
(Boehringer et al. 2001) or NORAS (Viana et al. 2001) , 
tracing the large scale structure (LSS) of the Universe up 
to <~ 1000 /i Mpc of distance (h is the Hubble constant in 
units of 100 km/s/Mpc). These data are consistent with 
the spectrum of a low-density CDM model, with matter 
density parameter Cl m ~ 0.3. Also an analysis of the red- 
shift dependence of the X-ray cluster luminosity function, 
based on a Press & Schechter approximation, sets a rela- 
tion between Q m and as, which requires again Q m < 0.5, 
0.6 (f2 m : matter density parameter; a$: squared mean 
density contrast, on the scale of 8/i _1 Mpc) (Borgani et 
al. 2001). Such constraints on the model will become 
more stringent, in the forthcoming years, thanks to XMM- 
Newton and Chandra data, providing high-z cluster sam- 
ples with tenths of thousands objects (see, e.g., Bartlett 
et al. 2001). 

In our opinion, such fast improvements of the observa- 
tional materials deserve to be matched by a suitable effort 
of the theoretical analysis. In this work we try to con- 
tribute to it, by providing an improved pattern to measure 
individual cluster densities, radii and masses in numerical 
simulations, without reference to a standard density con- 
trast. Although, both in real samples as well as in simu- 
lations, cluster members are directly observable together 
with their environs, a definition of cluster radii, masses, 
and density contrasts is needed both for statistical aims, 
when treating a large set of halos, and for individual ob- 
jects, e.g., to compare the observed mass distribution with 
analytical expressions of radial profiles. In the standard 
approach, a value for the density contrast is usually bor- 
rowed from analytical results, which are obtained under as- 
sumptions which can be incorrect in the real world, as well 
as in simulations. Such assumptions, however, were explic- 
itly considered to be the best possible choice. In particular, 
Cole & Lacey (1996) had explored the possibility to make 
use of the virial ratio, but considered its value in spheres, 
instead of layers. Here we see that, passing from the inte- 
gral to the differential virial ratio 72. (r) = —2T(r)/W(r), 
allows the detection of a neat boundary, where a transition 
from halo to surrounding material occurs. In fact, first of 
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all, (i) this inspection allows the detection of minimum 
points of 1Z in all halos, which could be easily missed in an 
integral inspection. Some of them are safely shown not to 
be numerical features, as they persist without appreciable 
shifts when we compare different resolution simulations of 
a given halo. Furthermore, (ii) we found that, CHb quite 
close to the radius f , where one of the persistent minima 
occurs, there is always a maximum of the function w(r); 
finally, CHb (iii) at the radius f , we have also that 

w(f) ~ K(f) . (6.1) 

Let us recall that, in general, w(r) is defined with refer- 
ence to the behavior of the potential energy around f, by 
setting: 

W(r) = W{f){r/f)~ w ^ (6.2) 

and, therefore, its r dependence accounts for the radial 
dependence of W(r). 

In this paper we have shown that the above properties 
(i), (ii), (iii) are consistent with the requirements that: (a) 
r sets a neat boundary between bound and unbound ma- 
terial, (b) A significant interval (depth ~ 50-100 /i _1 kpc), 
where w(r) is constant, lies about f. We called it a bound- 
ary layer (BL). Accordingly, the r dependence of W(r), 
inside of the BL, is expressed by eq. (6.1), by replacing 
the exponent w(r) with its constant value w. 

The relation between the observed features (i), (ii), (iii) 
and the properties (a), (b), was shown under the assump- 
tion of spherical symmetry. The fact that the conditions 
(i), (ii), (iii) are not satisfied exactly, but with a (slight) 
approximation, can be easily attributed to deviations of 
numerical structures from sphericity. However, such devi- 
ations ought to be quite mild, owing to the statistics we 
performed on a large number of halos, and this indicates 
that violent relaxation erased most previous a-spherical 
features. Non-spherical collapses, however, leave a sub- 
stantial imprint on the final density contrasts A c found for 
virialized halos whose spread accounts for different initial 
3-D geometries of fluctuations and different interactions 
with the environment during the collapsing stages. 

The simulations used to focus the above properties, at 
different levels of resolution, in TCDM and ACDM modes, 
were run both with the GADGET and the ART codes. 
Results do not depend on the code. We could also ex- 
tend the analysis, at the lowest resolution level, to sim- 
ulations run in quite large boxes (360 h^ 1 Mpc aside), 



where we had a large statistics of clusters (> 300 with 
M c > 4.2 • lO 14 ft. _1 M , for each model). We applied to 
all of them an i?w-requirement, meant to detect the BL, 
and found it in 97 % of the clusters considered. A de- 
tailed inspection of these clusters shows that the mass M c 
within the BL substantially coincides with the mass Md yn , 
evaluated from the velocities of all particles within r c — f, 
according to the virial theorem. Furthermore, a direct in- 
spection of the clusters without a BL, showed they were 
exceptional structures, mostly due to undergoing major 
merging, which certainly justified a dynamically unsettled 
condition and a badly a-spherical geometry. Let us how- 
ever notice (see Fig. 5) that "standard" a-sphericity of 
halo structures does not appear to be a serious obstacle to 
BL detection. 

At present, model clusters to fit X-ray data are mostly 
built using hydro-codes. This allows substantial improve- 
ments with respect to initial analyses based on pure New- 
tonian dynamics. There are, however, some features in 
the interpretation of both kinds of simulations, which rely 
on a sound definition of the virialization radius R c . This 
radius, whose setting immediately follows from the detec- 
tion of the BL position, corresponds to a wide spread of 
density contrasts A c , and, in general, is significantly dif- 
ferent from the radius R v , obtained using a given density 
contrast A„. Even keeping to models without hydrody- 
namics, it is clear that using different radii R c and density 
contrast A c could lead to different relations between Lx, 
Tx and M c . As a consequence the one-to-one correspon- 
dence between M c and Tx could be lost and it would no 
longer be granted that a mass function and a temperature 
function are substantially equivalent quantities. 

ACKNOWLEDGEMENTS 

We thank INAF for allowing us the CPU time to per- 
form the ART simulation C and D at the CINECA con- 
sortium (grant cnami44a on the SGI Origin 3800 ma- 
chine). GADGET simulations of high-resolution clusters 
have been run on the 16 Linux PC Beowulf cluster at the 
Osservatorio Astronomico di Torino. Alessandro Gardini 
is thanked for allowing us to make use of the simulations 
A and B; we also wish to thank him, Loris Colombo and 
Giuseppe Tormen for useful discussions during the prepa- 
ration of this paper. 



REFERENCES 



Bartlett J., et al., 2001 "Galaxy Clusters and the High Redshift 
Universe Observed in X-ravs" Proceedings of the XXI Moriond 
Conference, (preprint astro-ph/0106098 I 

Boehringer, H., et al., 2UUU, ApJS, izy, 435 B 

Borgani, S., et al., 2001, ApJ, 561, 13 B 

Brian, G. & Norman, M., 1998, ApJ, 495, 80 

Carlberg, R.G., et al., 1996, ApJ, 462, 32 

Cole, S. & Lacey, C, 1996, MNRAS, 281, 716 

Couchman, H. M. P., 1991, ApJ, 268, L23 

Doroshkevivh, A.G., Kotok, E., Novikov, I.,Polyudov, A.N., 

Shandarin, S.F. & Schctman, S., 1980, MNRAS, 192, 321 
Eke, V.R., Cole, S., & Frenk, C.S., 1996, MNRAS, 282, 263 
Eke, V., Navarro, J. & Frenk, C, 1998, ApJ, 503, 569E 
Gardini, A., Murante, G. & Bonometto, S. A., 1999, ApJ, 524, 510 
Ghigna S., Moore B., Governato F., Lake G., Quinn T. & Stadel J. 

, 2000, ApJ 544, 616 
Girardi M., Borgani S., Giuricin G., Mardirossian F. & Mezzetti M., 

1998, ApJ, 506, 45G 
Goldstein H., Poole C. & Safko, 2002, Classical Mechanics, 3rd ed. 

- Addison Wesley pub. 



Gott, R. & Rees, M. 1975, A&A, 45, 365G 

Governato F., Babul A., Quinn T., Tozzi P., Baugh CM., Katz N., 

& Lake G., 1999, MNRAS, 307, 949 
Kauffmann G., Colberg J., Diaferio A. & White S.D.M. 1999, 

MNRAS, 303, 188 
Klypin A., Kravtsov A., Bullock J. & Primack J., 2001. ApJ 554, 

903 

Knebe A., Kravtsov A., Gottlober S. & Klypin A., 2000 MNRAS, 
317, 630 

Kravtsov A., Klypin A. & Khokhlov A., 1997 ApJ, 111, 73 K 
Lacey, C.& Cole, S., 1993, MNRAS, 262, 627 
Lacey, C.& Cole, S., 1994, MNRAS, 271, 676 

Lahav, O., Lilje, P.R., Primack, J.R. & Rees, M., 1991, MNRAS, 
282, 263E 

Maccio A. V., Gardini A., Ghigna S. & Bonometto S. A., 2002, ApJ, 
564, 1M 

Monaco, P & Murante, G., Physical Review D , 60/10 , id. 103502 , 
1999 

Navarro J., Frenk C. & White S.D.M. 1996, ApJ, 462, 563 
Navarro J., Frenk C. & White S.D.M. 1997, ApJ, 490, 493 



10 



Mass of clusters in simulations 



Peebles P.J.E., 1980, The Large Scale Structure of the Universe, 

Princeton University Press, Princeton 
Power C, Navarro J., Jenkins A., Frenk C, White S.D.M., 

Springel V., Stadel J. & Quinn T., 2002, MNRAS (submitted) 

astro-ph/0201544 
Press w.ri. & scnecnter P., 1974, ApJ, 187, 425 
Sheth R.K., Mo H.J. & Tormen G., 2001 MNRAS, 323 ,1 



Sheth R.K. & Tormen G., 1999 MNRAS, 308, 119 

Sheth R.K. & Tormen G., 2002 MNRAS 329, 61 

Schueker P., Caldwell R.R., Boehringer H., Collins C.A. & Guzzo L., 

astro-ph/0211480 and A&A (in press) 
Springei v., Yosmaa N. & White S.D.M., 2001 New A, 6, 79 S 
Viana P., Nichol R. & Liddle A., 2002, Apj, 569L, 75, V 
Zel'dovich, Ya. B., 1970 A&A, 5,84 



Maccio et al. 



11 



APPENDIX A 

Let us assume that halos are spherically symmetric. If BL's are in virial equilibrium, at any r, in their interior 

(•>>--♦ - ii£l>, Ml) 

M(< r) being the mass inside r. According to kinetic gas theory, however, p = p(v 2 )/3. Accordingly, local equilibrium 
prescribes that p = —p4>/3 and, up to first order in Ar/r, 



3[p+U+ -p_U_] = —Gp_r 3 _-P[- 

r + 



Air 3-/3 r M+ M_ , 
a 7* . f 



3 2/3 - 1 M_ Ar 



~47rGp_ri— (A2) 

3 r_ r_ 

with a /3 value such that pr 3 ^ 13 is constant across Ar (indices ± indicate that a quantity is evaluated in r±). 
Let us compare this term with the potential energy of the layer 

W = G f + d»r,(r)^ * 4*Gp_r" ^ Ar 
j r _ r T— T— 



also evaluated assuming spherical symmetry and neglecting terms of higher order in Ar/r. The ratio between pV terms 
and potential energy is ~ (2/3 — l)/3 and would vanish for [3 — 0.5. 

Let us now assume to be dealing with halos whose profile is p(r) ~ p x( r / r s) with 

X (x) = [x a (l+xf- a }- 1 , (A4) 

where r s is the scale radius of the profile and the concentration of the halo is c = r s /r v . This profile is found for most 
halos in our simulations, up to r < r v , with a value of c in the interval 4-7 and values of a ~ 1; this agrees with previous 
numerical analyses (see, e.g., Navarro et al. 1996, 1997, Ghigna et al 2000, Klypin et al 2001, Power at al 2002), also 
extended to smaller values of r, finding values of a ranging between 0.8 and 1.2. Using this profile we can evaluate the 
r dependence of /3 for a values in the above interval. Such values are shown in Fig. 15 and indicate that for c ~ 4-7, 
(3 ~ 0.5-0.3. 

This shows that the contribution of the pressure terms to the virial balance is expected to be significantly smaller than 
the contribution of the potential term. 

Let us however add some further comments: 

(i) Using the profile in eq. (A4) for r >^ r v is an approximation. An inspection of numerical halos, already significantly 
inhomogeneous at such radii (see, e.g., Fig. 5), shows that the slope of the profile is often smoother there. 

(ii) Accordingly, owing to eqs. (A2)-(A3), the pV terms can be absorbed in a factor (~ 1), set in front of the potential 
term of the virial balance. This correction does not cause any displacement of the point where 1Z has its minima. On the 
contrary, it may interfere with the condition w = 1Z, which, in principle, can be suitably improved. 

(iii) Deviations from spherical symmetry, however, are far from being fully negligible and the discrepancies between 1Z 
and w, found in halo analysis, show that a correction of 1Z by ~ 10 %, would not improve the efficiency of the method. 

APPENDIX B 

We tried to find a BL, i.e. to apply the Rw requirement, also in two N-body simulations aimed at studying the evolution 
of a galactic stellar disk in a cosmological context. These simulations, performed with GADGET, will be discussed in a 
forthcoming paper (Murante, Curir & Mazzei, in preparation). Here we used simulations with DM only (no galactic disk 
is present). 

In the first one, initial conditions were set up with the same multi-mass technique described in Sec. 2, using the ART 
package. A DM halo, whose mass is M v = 1.14 x 10 11 /i" 1 M Q at redshift z = 0, was simulated at high-resolution using 
DM particles of mass M p = 1.21 x 10 6 /i _1 M Q . According to eq. (1.1), the virial radius r v w 0.125/i _1 Mpc. External 
forces were taken into account by setting heavier and heavier particles in 3 concentric shells. We checked that no intruders 
(heavier mass particles) were ever closer than 0.5 ft," 1 Mpc from the center of mass of the halo. The Rw requirement was 
applied to this halo, at redshifts z = and z = 2. 

At z — a BL was found, with a radius r c = 0.113/j" 1 Mpc, yielding a ratio M c /Md yn = 1.0502. On the contrary, at 
z=2, we could not find a good fit for the \ 2 parameter used to implement the Rw requirement. It should be noticed that 
these results held, in spite of the fact that the halo concentration is ~ 20. 

This suggests that BL's are found by the Rw requirement also when the condition pr 3 = const, is (mildly) violated; 
on the contrary, no BL is found when the the halo is not yet virialized, as it was at z=2 and is also confirmed by direct 
inspection. Apparently, the presence of a BL however distinguishes gravitationally relaxed objects from those still in an 
evolutionary phase. 
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A partial confirm of these two results comes from the application of the Rw requirement to a second simulation, following 
the evolution of an isolated halo of nominal mass equal to the previous one, with a nominal radius of R = 0.180 Mpc, 
generated with a NFW density profile (see eq. A4) and with initial particle velocities picked up from a multivariate 
Maxwellian distribution; the NFW density profile was cut off at the nominal radius; the concentration c was set to the 
value measured in the previous case (using r v to evaluate it, as in Navarro et al, 1997). The halo was evolved for w 10 Gyr, 
to reach the age of the Universe at z = 2 in a ACDM cosmology where Cl\ = 0.7 and Vt m = 0.3. 

We applied the Rw algorithm, to find a BL at the beginning (t = Gyr) and at the end (t = 10.23 Gyr) of the simulation. 
No BL could be found at t = 0, when the halo was not in equilibrium. Apparently, therefore, the presence of a BL is 
characteristic of gravitationally stable objects, independently of the geometrical distribution of particles. At t = 0, a SO 
algorithm, based on a purely geometrical prescription, would comfortably detect a halo. The Rw requirement, instead, 
which is based on dynamical prescriptions, appears to be more selective than purely geometrical recipes. 

At t = 10.23 Gyr, we find a layer fulfilling the Rw requirement, with r c = 0.023 Mpc, yielding M c /M dyn = 1.0201. This 
happens in spite of the value of the density contrast at r c , which is quite large (p/p cr it — 2457). In this (quite peculiar) 
case, the definition of virialization based on the density contrast fails, in spite of the fact that gravitational equilibrium is 
established. 

These arguments seem to indicate that a technique based on the Rw requirement is applicable on a wide range of DM 
halo masses, even when c is fairly large; mild deviations from the requirement that pr 3 is constant do not seem to destroy 
its effectiveness. Further analysis will be performed to put precise limits to the cases when the Rw criterion is a useful 
tool for telling apart virialized halos from still strongly evolving ones. 
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Parameters for the simulations performed by the AP3M code. 
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Fig. 1. — Integral (dashed line) and differential (solid line) virial ratios 1Z(> r) and TZ(r) are plotted against the halo radius r for a ACDM 
halo. Heavy dots indicate the points where the differential 1Z has minima. 
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Fig. 2. — The virial ratio 1Z = —2T/W is plotted against the halo radius r, at different resolution levels. The value of N shown in each 
box yields the number N 3 of particles used at the highest resolution level. Heavy dots show the setting of the permanent minima. The 
plots (a) and (b) refer to ACDM and TCDM clusters, respectively. The minima fulfilling the Rw requirement lay at ~ 1.62 h~ -'^Mpc and at 
~ 1.75 h" 1 Mpc for the halos a and b, respectively. 
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Fig. 3. — The TZ and w dependence on r is shown in the best and worst cases treated at the highest resolution level, (a and b, respectively); 
still at this resolution, the w behavior is rather noisy (compare however with Fig. 4, to see how better resolution reduces numerical noise). 
According to the Rw requirement, the case a shows a neat intersection of the two curves for a maximum of w and a minimum of 7Z. However, 
also in the case b, the correspondence between the maximum of w and the minimum of TZ is apparent. 
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Fig. 4. — A typical example of 1Z. and ui behaviors at the initial resolution level. This halo is the same shown in Fig. 3a, at higher resolution. 
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Fig. 5. — Image of a typical cluster, with the location of the boundary layer (BL). The cluster shown has a radius r c = 1.76 h~ 1 M.pc, a 
mass M c = 8.92 ■ 10 14 /i~ 1 AfQ and a density contrast A c = 141. The m.s. velocity of particles within r c is 1574km/s. This cluster is taken 
from simulation D. Gray scales refer to particle velocities. 





ACDM 

— A = 95.9 

,1 




i 1 


■ : i 



10 15 



Fig. 6. — The mass dependence of the density contrast is shown. Points refer to all clusters in the simulation B (ACDM). Heavy points 
give the average A c for each M c interval (constant logarithmic width). Bars are 1— <r errors of the averages. The thick horizontal line is the 
average value of A c . The dashed line is the virial density contrast A„ expected for a unperturbed spherical fluctuation growth. 
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Fig. 7. — The same as Fig. 6 for clusters in simulation A (TCDM). At high M c some heavy points are shown without error bars, when 
there is one cluster per logarithmic mass interval. 
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Fig. 8. — Average cluster radii (heavy dots) are plotted against their masses M c for the simulation A. The solid line is obtained for the 
average density contrast (150.5); the dashed line, instead, is obtained from the "virial" density contrast (178). Error bars yield the variance 
of cluster radii around their averages. 
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Fig. 9. — The same as Fig. 8, for the simulation B. Solid and dashed lines correspond to density contrasts 95.9 and 111, respectively. 
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Fig. 10. — We show how "virial" quantities can introduce spurious characteristic scales, if fit against "real" quantities. Here mean cluster 
radii (r c ) are plotted against M c (filled circles) and M v (empty circles). Error bars are omitted to avoid confusion (see, however, Fig. 7). The 
horizontal line could be interpreted as a mass— independent average radius for all clusters. 
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Fig. 11. — The ratio between masses evaluated from particle velocities, according to the virial theorem (M^jj), and summing particle 
masses (M c i uater ) is shown limiting clusters either at r c or at r v , for all clusters of simulation A. 
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Fig. 12. — The same as Fig. 11 for simulation B. 
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Fig. 15. — (i yields the deviation from 3 of the slope of the halo profile in the region where the BL is expected to lay. 



